#include "wrapper_hypre.h"

HYPRE_Solver par_solver, par_precond;
HYPRE_StructSolver s_solver, s_precond;
HYPRE_SStructSolver solver, precond;

void buildup_solver(std::string solver_name, std::string prec_name, const TEST_CONFIG & config_file)
{
    int my_pid; MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
    HYPRE_Int major, minor, patch, single;
    char * version;
    HYPRE_VersionNumber(&major, &minor, &patch, &single);
    HYPRE_Version(&version);
    if (my_pid == 0) printf("Major %d Minor %d Patch %d single %d :: %s\n", 
        major, minor, patch, single, version);
    // 设置求解器
    if      (prec_name == "AMG") {
        if (solver_name == "CG") {
            if (my_pid == 0) printf("using ParCSR PCG\n");
            HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &par_solver);
            HYPRE_ParCSRPCGSetMaxIter(par_solver, config_file.max_iter);
            HYPRE_ParCSRPCGSetTol(par_solver, config_file.rtol);
            HYPRE_ParCSRPCGSetAbsoluteTol(par_solver, config_file.atol);
            HYPRE_ParCSRPCGSetTwoNorm(par_solver, 1 );
            // HYPRE_ParCSRPCGSetRelChange(par_solver, 0 );
            HYPRE_ParCSRPCGSetPrintLevel(par_solver, config_file.iter_print_level);
            HYPRE_ParCSRPCGSetLogging(par_solver, 1);
        }
        else if (solver_name == "GMRES") {
            if (my_pid == 0) printf("using ParCSR GMRES\n");
            HYPRE_ParCSRGMRESCreate(MPI_COMM_WORLD, &par_solver);
            HYPRE_ParCSRGMRESSetKDim(par_solver, config_file.restart_len);
            HYPRE_ParCSRGMRESSetMaxIter(par_solver, config_file.max_iter);
            HYPRE_ParCSRGMRESSetTol(par_solver, config_file.rtol);
            HYPRE_ParCSRGMRESSetAbsoluteTol(par_solver, config_file.atol);
            
            HYPRE_ParCSRGMRESSetPrintLevel(par_solver, config_file.iter_print_level);
            HYPRE_ParCSRGMRESSetLogging(par_solver, 1);
        }
    }
    else if (prec_name == "PFMG" || prec_name == "SMG") {
        if (solver_name == "CG") {
            if (my_pid == 0) printf("using Struct PCG\n");
            HYPRE_StructPCGCreate(MPI_COMM_WORLD, &s_solver);
            HYPRE_StructPCGSetMaxIter(s_solver, config_file.max_iter);
            HYPRE_StructPCGSetTol(s_solver, config_file.rtol);
            HYPRE_StructPCGSetAbsoluteTol(s_solver, config_file.atol);
            HYPRE_StructPCGSetTwoNorm(s_solver, 1 );

            HYPRE_StructPCGSetPrintLevel(s_solver, config_file.iter_print_level);
            HYPRE_StructPCGSetLogging(s_solver, 1);
        }
        else if (solver_name == "GMRES") {
            if (my_pid == 0) printf("using Struct GMRES\n");
            HYPRE_StructGMRESCreate(MPI_COMM_WORLD, &s_solver);
            HYPRE_StructGMRESSetKDim(s_solver, config_file.restart_len);
            HYPRE_StructGMRESSetMaxIter(s_solver, config_file.max_iter);
            HYPRE_StructGMRESSetTol(s_solver, config_file.rtol);
            HYPRE_StructGMRESSetAbsoluteTol(s_solver, config_file.atol);
            
            HYPRE_StructGMRESSetPrintLevel(s_solver, config_file.iter_print_level);
            HYPRE_StructGMRESSetLogging(s_solver, 1);
        }
    }
    else if (solver_name == "CG") {
        if (my_pid == 0) printf("using SStruct PCG\n");
        HYPRE_SStructPCGCreate(MPI_COMM_WORLD, &solver);
        HYPRE_SStructPCGSetMaxIter(solver, config_file.max_iter);
        HYPRE_SStructPCGSetTol(solver, config_file.rtol);
        HYPRE_SStructPCGSetAbsoluteTol(solver, config_file.atol);
        HYPRE_SStructPCGSetTwoNorm(solver, 1 );
        // HYPRE_SStructPCGSetRelChange(solver, 0 );
        HYPRE_SStructPCGSetPrintLevel(solver, config_file.iter_print_level);
        HYPRE_SStructPCGSetLogging(solver, 1);
    }
    else if (solver_name == "GMRES") {
        if (my_pid == 0) printf("using SStruct GMRES\n");
        HYPRE_SStructGMRESCreate(MPI_COMM_WORLD, &solver);
        HYPRE_SStructGMRESSetKDim(solver, config_file.restart_len);
        HYPRE_SStructGMRESSetMaxIter(solver, config_file.max_iter);
        HYPRE_SStructGMRESSetTol(solver, config_file.rtol);
        HYPRE_SStructGMRESSetAbsoluteTol(solver, config_file.atol);
        
        HYPRE_SStructGMRESSetPrintLevel(solver, config_file.iter_print_level);
        HYPRE_SStructGMRESSetLogging(solver, 1);
    }
    else MPI_Abort(MPI_COMM_WORLD, -20240123);

    // 设置预条件子
    if      (strstr(prec_name.c_str(), "Split")) {
        if (my_pid == 0) printf("using SStruct Split\n");
        HYPRE_SStructSplitCreate(MPI_COMM_WORLD, &precond);
        HYPRE_SStructSplitSetMaxIter(precond, 1);
        HYPRE_SStructSplitSetTol(precond, 0.0);
        HYPRE_SStructSplitSetPrintLevel(precond, config_file.prec_print_level);
        HYPRE_SStructSplitSetLogging(precond, 0);

        HYPRE_SStructSplitSetZeroGuess(precond);
        if      (prec_name == "Split-SMG" ) HYPRE_SStructSplitSetStructSolver(precond, HYPRE_SMG);
        else if (prec_name == "Split-PFMG") HYPRE_SStructSplitSetStructSolver(precond, HYPRE_PFMG);
        else MPI_Abort(MPI_COMM_WORLD, -1010);
        if (solver_name == "CG") 
            HYPRE_SStructPCGSetPrecond  (solver, HYPRE_SStructSplitSolve, HYPRE_SStructSplitSetup, precond);
        else if (solver_name == "GMRES")
            HYPRE_SStructGMRESSetPrecond(solver, HYPRE_SStructSplitSolve, HYPRE_SStructSplitSetup, precond);
    }
    else if (prec_name == "PFMG") {
        if (my_pid == 0) printf("using Struct PFMG\n");
        HYPRE_StructPFMGCreate(MPI_COMM_WORLD, &s_precond);
        HYPRE_StructPFMGSetMaxIter(s_precond, 1);
        HYPRE_StructPFMGSetTol(s_precond, 0.0);
        HYPRE_StructPFMGSetPrintLevel(s_precond, config_file.prec_print_level);
        HYPRE_StructPFMGSetLogging(s_precond, 0);
        
        HYPRE_StructPFMGSetZeroGuess(s_precond);
        if (solver_name == "CG") 
            HYPRE_StructPCGSetPrecond  (s_solver, HYPRE_StructPFMGSolve, HYPRE_StructPFMGSetup, s_precond);
        else if (solver_name == "GMRES")
            HYPRE_StructGMRESSetPrecond(s_solver, HYPRE_StructPFMGSolve, HYPRE_StructPFMGSetup, s_precond);
    }
    else if (prec_name == "SMG") {
        if (my_pid == 0) printf("using Struct SMG\n");
        HYPRE_StructSMGCreate(MPI_COMM_WORLD, &s_precond);
        HYPRE_StructSMGSetMaxIter(s_precond, 1);
        HYPRE_StructSMGSetTol(s_precond, 0.0);
        HYPRE_StructSMGSetPrintLevel(s_precond, config_file.prec_print_level);
        HYPRE_StructSMGSetLogging(s_precond, 0);
        
        HYPRE_StructSMGSetZeroGuess(s_precond);
        if (solver_name == "CG") 
            HYPRE_StructPCGSetPrecond  (s_solver, HYPRE_StructSMGSolve, HYPRE_StructSMGSetup, s_precond);
        else if (solver_name == "GMRES")
            HYPRE_StructGMRESSetPrecond(s_solver, HYPRE_StructSMGSolve, HYPRE_StructSMGSetup, s_precond);
    }
    else if (strstr(prec_name.c_str(), "SSAMG")) {
        if (my_pid == 0) printf("using SStruct SSAMG\n");
        HYPRE_SStructSSAMGCreate(MPI_COMM_WORLD, &precond);
        HYPRE_SStructSSAMGSetMaxIter(precond, 1);
        HYPRE_SStructSSAMGSetTol(precond, 0.0);
        HYPRE_SStructSSAMGSetPrintLevel(precond, config_file.prec_print_level);
        HYPRE_SStructSSAMGSetLogging(precond, 0);

        HYPRE_SStructSSAMGSetZeroGuess(precond);
        // 0: Jacobi, 1: Weighted Jacobi (default), 2: L1-Jacobi,
        // 10: Red/Black Gauss-Seidel (symmetric: RB pre-relaxation, BR post-relaxation)
        // HYPRE_SStructSSAMGSetRelaxType(precond, 10);// 10的话会报错

        if (prec_name == "SSAMG-base") {
            // HYPRE_SStructSSAMGSetMaxLevels(precond, 25);// 半结构化部分的层数，最粗层可以转到AMG或者weighted Jacobi
            // HYPRE_SStructSSAMGSetCoarseSolverType(precond, 1);// 0: weighted Jacobi, 1: BoomerAMG
        } else if (prec_name == "SSAMG-skip") {
            // HYPRE_SStructSSAMGSetMaxLevels(precond, 25);// 半结构化部分的层数，最粗层可以转到AMG或者weighted Jacobi
            // HYPRE_SStructSSAMGSetCoarseSolverType(precond, 1);// 0: weighted Jacobi, 1: BoomerAMG
            HYPRE_SStructSSAMGSetSkipRelax(precond, 1);
        } else if (prec_name == "SSAMG-hybrid") {
            HYPRE_SStructSSAMGSetMaxLevels(precond, 10);// 半结构化部分的层数，最粗层可以转到AMG
            HYPRE_SStructSSAMGSetCoarseSolverType(precond, 1);// 0: weighted Jacobi, 1: BoomerAMG
        } else if (prec_name == "SSAMG-opt") {
            HYPRE_SStructSSAMGSetMaxLevels(precond, 7);// 半结构化部分的层数，最粗层可以转到AMG
            HYPRE_SStructSSAMGSetCoarseSolverType(precond, 1);// 0: weighted Jacobi, 1: BoomerAMG
        } else MPI_Abort(MPI_COMM_WORLD, -20240125);

        if (config_file.case_name == "MBSM_bench" || config_file.case_name == "SAMR_bench") {
            HYPRE_SStructSSAMGSetRelaxWeight(precond, 0.6667);// weighted L1-Jacobi smoother factor equal to 3/2. ??? 迭代法运行2步就停止了
        }
        // 天气预报算例
        else if (config_file.case_name == "gfs") {
            // HYPRE_SStructSSAMGSetNumPreRelax(precond, 2);
            // HYPRE_SStructSSAMGSetNumPostRelax(precond, 2);
            HYPRE_SStructSSAMGSetMaxLevels(precond, 4);

            // HYPRE_SStructSSAMGSetRelaxType(precond, 0);// 0 和 1好像没区别，2会更慢一些，10也很慢
            // HYPRE_SStructSSAMGSetNonGalerkinRAP(precond, 1);// 如果设置为非galerkin会报错！

            // HYPRE_SStructSSAMGSetMaxCoarseSize(precond, 40000000);// 3层网格，56次迭代 （This is BEST）
            // HYPRE_SStructSSAMGSetMaxCoarseSize(precond, 20000000);// 4层网格，收敛明显变慢，即使光滑权重仍一直为0.67
            // HYPRE_SStructSSAMGSetMaxCoarseSize(precond, 10000000);// 5层网格
            // HYPRE_SStructSSAMGSetRelaxWeight(precond, 0.8);

            // 所建立的多层网格的信息跟通过SetMaxCoarseSize(40000000)的完全一致，但是最粗层的光滑权重却不一样，
            // 导致收敛明显变慢
            // HYPRE_SStructSSAMGSetMaxLevels(precond, 3);

            HYPRE_SStructSSAMGSetNumCoarseRelax(precond, 1);// 设置最粗层只做1次平滑
            // HYPRE_SStructSSAMGSetCoarseSolverType(precond, 1);// 即使在最好的SSAMG配置（3层，权重0.67）下将最粗层转为BoomerAMG也会完全不收敛
        }
        // 轮船制造算例
        else if (config_file.case_name == "kvlcc2") {
            
        }
        // 油藏模拟算例
        else if (config_file.case_name == "oil") {
        }
        // 海洋模拟算例
        else if (config_file.case_name == "ocean") {
            assert(false);
            // not used
        }

        // HYPRE_SStructSSAMGSetNonGalerkinRAP(precond, 1);// 会报错


        if (solver_name == "CG") 
            HYPRE_SStructPCGSetPrecond  (solver, HYPRE_SStructSSAMGSolve, HYPRE_SStructSSAMGSetup, precond);
        else if (solver_name == "GMRES")
            HYPRE_SStructGMRESSetPrecond(solver, HYPRE_SStructSSAMGSolve, HYPRE_SStructSSAMGSetup, precond);
    }
    else if (prec_name == "AMG") {
        if (my_pid == 0) printf("using ParCSR BoomerAMG\n");
        HYPRE_BoomerAMGCreate(&par_precond);
        HYPRE_BoomerAMGSetMaxIter(par_precond, 1);
        HYPRE_BoomerAMGSetTol(par_precond, 0.0);
        HYPRE_BoomerAMGSetPrintLevel(par_precond, config_file.prec_print_level);
        HYPRE_BoomerAMGSetLogging(par_precond, 0);        
        HYPRE_BoomerAMGSetMaxLevels(par_precond, 25);// 对于kvlcc2当变成5层时才有突变

        HYPRE_BoomerAMGSetCoarsenType(par_precond, 10);
        HYPRE_BoomerAMGSetAggNumLevels(par_precond, 1);//  one level (first) of aggressive coarsening
        HYPRE_BoomerAMGSetAggInterpType(par_precond, 4);// with multi-pass interpolation

        // 理想算例
        if (config_file.case_name == "MBSM_bench" || config_file.case_name == "SAMR_bench") {
            HYPRE_BoomerAMGSetRelaxType(par_precond, 8);// Forward/Backward L1-Gauss-Seidel 
            HYPRE_BoomerAMGSetStrongThreshold(par_precond, 0.25);

            // HYPRE_BoomerAMGSetInterpType(par_precond, 17);// in the following levels, matrix-based extended+i interpolation [23]
            HYPRE_BoomerAMGSetPMaxElmts(par_precond, 4);// truncated to a maximum of four nonzero coefficients per row.
        }
        // 天气预报算例
        else if (config_file.case_name == "gfs") {
            HYPRE_BoomerAMGSetSmoothType(par_precond, 5);// 用ILU作为复杂平滑子
            HYPRE_BoomerAMGSetSmoothNumLevels(par_precond, 5);// 前多少层都用ILU作为平滑子

            HYPRE_BoomerAMGSetRelaxType(par_precond, 6);// Forward/Backward L1-Gauss-Seidel （剩余层的平滑子）
            HYPRE_BoomerAMGSetStrongThreshold(par_precond, 0.25);
        }
        // 轮船制造算例
        else if (config_file.case_name == "kvlcc2") {
            
        }
        // 油藏模拟算例
        else if (config_file.case_name == "oil") {
            // HYPRE_BoomerAMGSetSmoothType(par_precond, 5);
            // HYPRE_BoomerAMGSetSmoothNumLevels(par_precond, 4);

            // HYPRE_BoomerAMGSetRelaxType(par_precond, 6);// Forward/Backward L1-Gauss-Seidel
            HYPRE_BoomerAMGSetStrongThreshold(par_precond, 0.25);
            // HYPRE_BoomerAMGSetPMaxElmts(par_precond, 5);
            // HYPRE_BoomerAMGSetRelaxOrder(par_precond, 1);// CF order according to csr.fasp
        }
        // 海洋模拟算例
        else if (config_file.case_name == "ocean") {
            HYPRE_BoomerAMGSetRelaxType(par_precond, 6);// Forward/Backward L1-Gauss-Seidel
            HYPRE_BoomerAMGSetStrongThreshold(par_precond, 0.25);
        }
        HYPRE_BoomerAMGSetModuleRAP2(par_precond, 1);// modularized option for computing the Galerkin product RAP

        if (solver_name == "CG") 
            HYPRE_ParCSRPCGSetPrecond  (par_solver, HYPRE_BoomerAMGSolve, HYPRE_BoomerAMGSetup, par_precond);
        else if (solver_name == "GMRES")
            HYPRE_ParCSRGMRESSetPrecond(par_solver, HYPRE_BoomerAMGSolve, HYPRE_BoomerAMGSetup, par_precond);
    }
    else {
        if (my_pid == 0) printf("No Precond!\n");
    }
}

void setup_and_solve(std::string solver_name, std::string prec_name, void* A, void* b, void* x,
    HYPRE_Real & t_setup, HYPRE_Real & t_solve, HYPRE_Real & final_res_norm, HYPRE_Int & num_iterations)
{
    // 开始求解
    if      (prec_name == "AMG") {
        if (solver_name == "CG") {
            MPI_Barrier(MPI_COMM_WORLD);
            t_setup = MPI_Wtime();
            HYPRE_ParCSRPCGSetup(par_solver, (HYPRE_ParCSRMatrix) A, (HYPRE_ParVector) b, (HYPRE_ParVector) x );
            t_setup = MPI_Wtime() - t_setup;
            
            MPI_Barrier(MPI_COMM_WORLD);
            t_solve = MPI_Wtime();
            HYPRE_ParCSRPCGSolve(par_solver, (HYPRE_ParCSRMatrix) A, (HYPRE_ParVector) b, (HYPRE_ParVector) x );
            t_solve = MPI_Wtime() - t_solve;
            
            HYPRE_ParCSRPCGGetNumIterations( par_solver, &num_iterations );
            HYPRE_ParCSRPCGGetFinalRelativeResidualNorm( par_solver, &final_res_norm );
        }
        else if (solver_name == "GMRES") {
            MPI_Barrier(MPI_COMM_WORLD);
            t_setup = MPI_Wtime();
            HYPRE_ParCSRGMRESSetup(par_solver, (HYPRE_ParCSRMatrix) A, (HYPRE_ParVector) b, (HYPRE_ParVector) x );
            t_setup = MPI_Wtime() - t_setup;
            
            MPI_Barrier(MPI_COMM_WORLD);
            t_solve  = MPI_Wtime();
            HYPRE_ParCSRGMRESSolve(par_solver, (HYPRE_ParCSRMatrix) A, (HYPRE_ParVector) b, (HYPRE_ParVector) x );
            t_solve = MPI_Wtime() - t_solve;
            HYPRE_ParCSRGMRESGetNumIterations( par_solver, &num_iterations );
            HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm( par_solver, &final_res_norm );
        }
    }
    
    else if (prec_name == "PFMG" || prec_name == "SMG") {
        if (solver_name == "CG") {
            MPI_Barrier(MPI_COMM_WORLD);
            t_setup = MPI_Wtime();
            HYPRE_StructPCGSetup(s_solver, (HYPRE_StructMatrix) A, (HYPRE_StructVector) b, (HYPRE_StructVector) x );
            t_setup = MPI_Wtime() - t_setup;
            
            MPI_Barrier(MPI_COMM_WORLD);
            t_solve = MPI_Wtime();
            HYPRE_StructPCGSolve(s_solver, (HYPRE_StructMatrix) A, (HYPRE_StructVector) b, (HYPRE_StructVector) x );
            t_solve = MPI_Wtime() - t_solve;
            
            HYPRE_StructPCGGetNumIterations( s_solver, &num_iterations );
            HYPRE_StructPCGGetFinalRelativeResidualNorm( s_solver, &final_res_norm );
        }
        else if (solver_name == "GMRES") {
            MPI_Barrier(MPI_COMM_WORLD);
            t_setup = MPI_Wtime();
            HYPRE_StructGMRESSetup(s_solver, (HYPRE_StructMatrix) A, (HYPRE_StructVector) b, (HYPRE_StructVector) x );
            t_setup = MPI_Wtime() - t_setup;
            
            MPI_Barrier(MPI_COMM_WORLD);
            t_solve  = MPI_Wtime();
            HYPRE_StructGMRESSolve(s_solver, (HYPRE_StructMatrix) A, (HYPRE_StructVector) b, (HYPRE_StructVector) x );
            t_solve = MPI_Wtime() - t_solve;
            HYPRE_StructGMRESGetNumIterations( s_solver, &num_iterations );
            HYPRE_StructGMRESGetFinalRelativeResidualNorm( s_solver, &final_res_norm );
        }
    }
    else if (solver_name == "CG") {
        MPI_Barrier(MPI_COMM_WORLD);
        t_setup = MPI_Wtime();
        HYPRE_SStructPCGSetup(solver, (HYPRE_SStructMatrix) A, (HYPRE_SStructVector) b, (HYPRE_SStructVector) x );
        t_setup = MPI_Wtime() - t_setup;
        
        MPI_Barrier(MPI_COMM_WORLD);
        t_solve = MPI_Wtime();
        HYPRE_SStructPCGSolve(solver, (HYPRE_SStructMatrix) A, (HYPRE_SStructVector) b, (HYPRE_SStructVector) x );
        t_solve = MPI_Wtime() - t_solve;
        
        HYPRE_SStructPCGGetNumIterations( solver, &num_iterations );
        HYPRE_SStructPCGGetFinalRelativeResidualNorm( solver, &final_res_norm );
    }
    else if (solver_name == "GMRES") {
        MPI_Barrier(MPI_COMM_WORLD);
        t_setup = MPI_Wtime();
        HYPRE_SStructGMRESSetup(solver, (HYPRE_SStructMatrix) A, (HYPRE_SStructVector) b, (HYPRE_SStructVector) x );
        t_setup = MPI_Wtime() - t_setup;

        MPI_Barrier(MPI_COMM_WORLD);
        t_solve  = MPI_Wtime();
        HYPRE_SStructGMRESSolve(solver, (HYPRE_SStructMatrix) A, (HYPRE_SStructVector) b, (HYPRE_SStructVector) x );
        t_solve = MPI_Wtime() - t_solve;

        HYPRE_SStructGMRESGetNumIterations( solver, &num_iterations );
        HYPRE_SStructGMRESGetFinalRelativeResidualNorm( solver, &final_res_norm );
    }
}

void destroy_solver(std::string solver_name, std::string prec_name)
{
    if (prec_name == "AMG") {
        // stat_amg_pattern(prec_name);
        HYPRE_BoomerAMGDestroy(par_precond);
        if      (solver_name == "CG")    HYPRE_ParCSRPCGDestroy(par_solver);
        else if (solver_name == "GMRES") HYPRE_ParCSRGMRESDestroy(par_solver);
    } 
    else if (prec_name == "PFMG" || prec_name == "SMG") {
        if (prec_name == "PFMG") HYPRE_StructPFMGDestroy(s_precond);
        else                     HYPRE_StructSMGDestroy(s_precond);
        if      (solver_name == "CG")    HYPRE_StructPCGDestroy(s_solver);
        else if (solver_name == "GMRES") HYPRE_StructGMRESDestroy(s_solver);
    }
    else {
        if      (strstr(prec_name.c_str(), "Split")) HYPRE_SStructSplitDestroy(precond);
        else if (strstr(prec_name.c_str(), "SSAMG")) HYPRE_SStructSSAMGDestroy(precond);
        if      (solver_name == "CG")    HYPRE_SStructPCGDestroy(solver);
        else if (solver_name == "GMRES") HYPRE_SStructGMRESDestroy(solver);
    }
}

// HYPRE_Real check_residual(const HYPRE_ParCSRMatrix parcsr_A, const HYPRE_ParVector par_x,
//     const HYPRE_ParVector par_b, HYPRE_ParVector par_y)
// {
//     // 计算真实残差
//     HYPRE_ParVectorCopy(par_b, par_y);// y = b
//     HYPRE_ParCSRMatrixMatvec(-1.0, parcsr_A, par_x, 1.0, par_y);// y += -A*x
//     HYPRE_Real r_nrm2 = 0.0, b_nrm2 = 0.0;
//     HYPRE_ParVectorInnerProd(par_b, par_b, &b_nrm2); b_nrm2 = sqrt(b_nrm2);
//     HYPRE_ParVectorInnerProd(par_y, par_y, &r_nrm2); r_nrm2 = sqrt(r_nrm2);
//     int my_pid; MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
//     if (my_pid == 0) printf("||r||_2 = %.7e    ||b||_2 = %.7e    ||r||_2/||b||_2 = %.7e\n", r_nrm2, b_nrm2, r_nrm2/b_nrm2);
//     return r_nrm2;
// }

void check_residual(std::string prec_name,
    const HYPRE_SStructMatrix ss_A, const HYPRE_SStructVector ss_x,
    const HYPRE_SStructVector ss_b, HYPRE_SStructVector ss_y, HYPRE_Real & r_nrm2, HYPRE_Real & b_nrm2)
{
    // 计算真实残差
    if (prec_name == "PFMG" || prec_name == "SMG") {
        HYPRE_StructMatrix s_A;
        HYPRE_StructVector s_x, s_b, s_y;
        HYPRE_SStructMatrixGetObject(ss_A, (void **) &s_A);
        HYPRE_SStructVectorGetObject(ss_b, (void **) &s_b);
        HYPRE_SStructVectorGetObject(ss_x, (void **) &s_x);
        HYPRE_SStructVectorGetObject(ss_y, (void **) &s_y);

        HYPRE_StructVectorCopy( s_b, s_y);
        HYPRE_StructMatrixMatvec( -1.0, s_A, s_x, 1.0, s_y);
        HYPRE_StructVectorInnerProd( s_b, s_b, &b_nrm2); b_nrm2 = sqrt(b_nrm2);
        HYPRE_StructVectorInnerProd( s_y, s_y, &r_nrm2); r_nrm2 = sqrt(r_nrm2);
    } else {
        HYPRE_SStructVectorCopy(ss_b, ss_y);// y = b
        HYPRE_SStructMatrixMatvec(-1.0, ss_A, ss_x, 1.0, ss_y);// y += -A*x
        HYPRE_SStructInnerProd(ss_b, ss_b, &b_nrm2); b_nrm2 = sqrt(b_nrm2);
        HYPRE_SStructInnerProd(ss_y, ss_y, &r_nrm2); r_nrm2 = sqrt(r_nrm2);
    }
}

// void check_residual(const HYPRE_SStructMatrix ss_A, const HYPRE_SStructVector ss_x,
//     const HYPRE_SStructVector ss_b, HYPRE_SStructVector ss_y, HYPRE_Real & r_nrm2, HYPRE_Real & b_nrm2)
// {
//     // 计算真实残差
//     HYPRE_SStructVectorCopy(ss_b, ss_y);// y = b
//     HYPRE_SStructMatrixMatvec(-1.0, ss_A, ss_x, 1.0, ss_y);// y += -A*x
//     HYPRE_SStructInnerProd(ss_b, ss_b, &b_nrm2); b_nrm2 = sqrt(b_nrm2);
//     HYPRE_SStructInnerProd(ss_y, ss_y, &r_nrm2); r_nrm2 = sqrt(r_nrm2);
//     // int my_pid; MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
//     // if (my_pid == 0) printf("||r||_2 = %.7e    ||b||_2 = %.7e    ||r||_2/||b||_2 = %.7e\n", r_nrm2, b_nrm2, r_nrm2/b_nrm2);
// }

void stat_amg_pattern(std::string prec_name)
{
    if (prec_name != "AMG") return ;

    // int glb_pid; MPI_Comm_rank(MPI_COMM_WORLD, & glb_pid);

    hypre_ParAMGData * amg_data = (hypre_ParAMGData*) par_precond;

    
    const HYPRE_Int num_levels = amg_data->num_levels;
    // if (glb_pid == 0) printf("Num levels %d\n", num_levels);
    for (int ilev = 0; ilev < num_levels; ilev ++) {
        const hypre_ParCSRMatrix * par_A = amg_data->A_array[ilev];
        // if (glb_pid == 0) printf("%d nrows %d\n", ilev, par_A->global_num_rows);

        const hypre_ParCSRCommPkg * comm_pkg = par_A->comm_pkg;
        // printf("lev %d Proc %d num_recv %d num_send %d\n", ilev, glb_pid, comm_pkg->num_recvs, comm_pkg->num_sends);

        const MPI_Comm lev_comm = comm_pkg->comm;
        int lev_pid, lev_nprocs;
        MPI_Comm_rank(lev_comm, & lev_pid);
        MPI_Comm_size(lev_comm, & lev_nprocs);

        // if (glb_pid == 0) printf("lev %d num_Procs %d\n", ilev, lev_nprocs);

        FILE * fp = nullptr;
        if (lev_pid == 0) { fp = fopen(("amg.L"+std::to_string(ilev)).c_str(), "w"); fclose(fp); }

        for (int p = 0; p < lev_nprocs; p++) {
            if (p == lev_pid) {
                fp = fopen(("amg.L"+std::to_string(ilev)).c_str(), "a");// attach

                // for (int is = 0; is < comm_pkg->num_sends; is++) {
                //     const int dst_pid = comm_pkg->send_procs[is];
                //     const int num = comm_pkg->send_map_starts[is + 1] - comm_pkg->send_map_starts[is];
                //     printf("lev %d Proc %d send  to  %d num %d\n", ilev, glb_pid, dst_pid, num);
                // }

                for (int ir = 0; ir < comm_pkg->num_recvs; ir++) {
                    const int src_pid = comm_pkg->recv_procs[ir];
                    const int num = comm_pkg->recv_vec_starts[ir + 1] - comm_pkg->recv_vec_starts[ir];
                    // printf("lev %d Proc %d recv from %d num %d\n", ilev, glb_pid, src_pid, num);
                    fprintf(fp, "%d %d %d\n", lev_pid, src_pid, num);
                }
                fclose(fp);
                fflush(stdout);
            }
            MPI_Barrier(lev_comm);
        }// process print sequentially
    }// lev loop
}